#
# 主逻辑代码
# 

# 加载代码所需的环境依赖库
import os
import sys
import logging
from logHandler import LogHandler
LogHandler.init_log_handler('run_log\\'+"leafArea")
logger = logging.getLogger("mylog")
if os.system('python setup.py build_ext --inplace')!=0: # 在当前环境下执行此命令
    logger.info(" 没有正确编译pyx 文件，退出程序")
    sys.exit(-1)

import math
import numpy as np
import random
import scipy
# 最小二乘求解非线性方程组
from scipy.optimize import leastsq,fsolve,root
from osgeo import gdal,gdalconst
import pandas as pds
from scipy import interpolate
from multiprocessing import pool

# 加载自己库
from ImageHandle import ImageHandler
from sample_process import read_sample_csv,combine_sample_attr,ReprojectImages2,read_tiff,check_sample,split_sample_list
from LAIProcess import train_WMCmodel,test_WMCModel,process_tiff

from model_evaluate import model_evaluate

# 常量声明区域
imageHandler=ImageHandler()


if __name__=="__main__":
    # 输入
    # a, 工作空间
    work_path = r"D:\BaiduNetdiskDownload\MicroWorkspace\data\temp"
    # b. 结果工作
    result_dir_path = r"D:\BaiduNetdiskDownload\MicroWorkspace\data\result"
    # c. 是否绘制图像
    draw_flag = True
    # 1. 后向散射系数 dB
    sigma_path = r"D:\BaiduNetdiskDownload\MicroWorkspace\data\锡林浩特数据\hv_pwr_fill_db.tif"
    # 2. 局地入射角
    incident_angle_path = r"D:\BaiduNetdiskDownload\MicroWorkspace\data\LocalIncidenceAngle.tif"
    # 3. 样本csv地址
    lai_csv_path = r"D:\BaiduNetdiskDownload\MicroWorkspace\laiplus2.csv"
    # 4. NDVI影像地址 -- 修正模型
    NDVI_tiff_path = r"D:\BaiduNetdiskDownload\MicroWorkspace\data\锡林浩特数据\202208NDVI.tif"
    # 5. 土壤含水量影像地址
    soil_water_tiff_path = r"D:\BaiduNetdiskDownload\MicroWorkspace\data\202208.tif"
    
    soil_simga_path=r"D:\BaiduNetdiskDownload\MicroWorkspace\data\dbl_grid.tif"
    # 6. 土壤含水量样本地址
    soil_water_csv_path = r""

    # 7. 选择土壤含水量影像
    soil_water = 'tiff'
    # 8. 输出图片
    train_err_image_path=os.path.join(result_dir_path,"train_image.png")
    test_err_image_path=os.path.join(result_dir_path,"test_image.png")
    evalution_image_path=os.path.join(result_dir_path,"evaluation_image.png")

    # 临时变量
    soil_tiff_resample_path=os.path.join(work_path,"soil_water.tiff")   # 与 后向散射系数同样分辨率的 土壤水分影像
    NDVI_tiff_resample_path=os.path.join(work_path,'NDVI.tiff') # 与 后向散射系数产品同样分辨率的 NDVI影像
    incident_angle_resample_path=os.path.join(work_path,"localincangle.tiff") 
    soil_simga_dbl_path=os.path.join(work_path,"dbl_temp.tiff")
    # 读取数据
    lai_sample=read_sample_csv(lai_csv_path)        # 读取样本数据
    sigma_tiff=read_tiff(sigma_path)                # 读取后向散射系数
    incident_angle=read_tiff(incident_angle_path)   # 读取局地入射角
    
    # 对于土壤水分、NDVI做重采样
    ReprojectImages2(soil_water_tiff_path,sigma_path,soil_tiff_resample_path,resampleAlg=gdalconst.GRA_Bilinear) 
    ReprojectImages2(NDVI_tiff_path,sigma_path,NDVI_tiff_resample_path,resampleAlg=gdalconst.GRA_Bilinear) 
    ReprojectImages2(incident_angle_path,sigma_path,incident_angle_resample_path,resampleAlg=gdalconst.GRA_NearestNeighbour) 
    ReprojectImages2(soil_simga_path,sigma_path,soil_simga_dbl_path,resampleAlg=gdalconst.GRA_NearestNeighbour) 
    
    
    soil_water_tiff=read_tiff(soil_tiff_resample_path) # 读取土壤含水量影像
    NDVI_tiff=read_tiff(NDVI_tiff_resample_path) # 引入NDVI
    incident_angle=read_tiff(incident_angle_resample_path)   # 读取局地入射角
    soil_simga_dbl=read_tiff(soil_simga_dbl_path)
    
    soil_water_tiff['data']=soil_water_tiff['data']/100.0       # 转换为百分比
    incident_angle['data']=incident_angle['data']*np.pi/180.0   # 转换为弧度值
    sigma_tiff['data']=np.power(10,(sigma_tiff['data']/10))     # 转换为线性值
    
    # float32 转 float64
    soil_water_tiff['data']=soil_water_tiff['data'].astype(np.float64)
    incident_angle['data']=incident_angle['data'].astype(np.float64)
    sigma_tiff['data']=sigma_tiff['data'].astype(np.float64)    
    soil_simga_dbl['data']=soil_simga_dbl['data'].astype(np.float64)
    # 将土壤水分与lai样本之间进行关联
    lai_water_sample=[] # ['日期', '样方编号', '经度', '纬度', 'LAI','土壤含水量']
    if soil_water=='tiff': 
        lai_water_sample=combine_sample_attr(lai_sample,soil_water_tiff)
        pass
    else: # 这个暂时没有考虑
        pass
    
    # 将入射角、后向散射系数与lai样本之间进行关联
    lai_water_inc_list=combine_sample_attr(lai_water_sample,incident_angle)         # ['日期','样方编号','经度','纬度','叶面积指数',"后向散射系数",'土壤含水量','入射角']
    lai_waiter_inc_sigma_list=combine_sample_attr(lai_water_inc_list,sigma_tiff)    # ['日期','样方编号','经度','纬度','叶面积指数',"后向散射系数",'土壤含水量','入射角','后向散射系数']   
    lai_waiter_inc_sigma_list=combine_sample_attr(lai_waiter_inc_sigma_list,soil_simga_dbl)    # ['日期','样方编号','经度','纬度','叶面积指数',"后向散射系数",'土壤含水量','入射角','后向散射系数','soil_simga_dbl']    
    #lai_waiter_inc_sigma_NDVI_list=combine_sample_attr(lai_waiter_inc_sigma_list,NDVI_tiff)    # ['日期','样方编号','经度','纬度','叶面积指数',"后向散射系数",'土壤含水量','入射角','后向散射系数','NDVI']    
    lai_waiter_inc_sigma_list=check_sample(lai_waiter_inc_sigma_list)               # 清理样本 ['日期','样方编号','经度','纬度','叶面积指数',"后向散射系数",'土壤含水量','入射角','后向散射系数'] 
    #lai_waiter_inc_sigma_NDVI_list=check_sample(lai_waiter_inc_sigma_NDVI_list) # ['日期','样方编号','经度','纬度','叶面积指数',"后向散射系数",'土壤含水量','入射角','后向散射系数','NDVI'] 
    # 数据集筛选
    lai_waiter_inc_sigma_list_result=[]
    # 筛选保留的数据集
    logger.info("保留得数据集如下")
    for i in range(len(lai_waiter_inc_sigma_list)):
        if i in []:
            continue
        logger.info(str(lai_waiter_inc_sigma_list[i]))
        lai_waiter_inc_sigma_list_result.append(lai_waiter_inc_sigma_list[i])
    lai_waiter_inc_sigma_list=lai_waiter_inc_sigma_list_result

    #[sample_train,sample_test]=split_sample_list(lai_waiter_inc_sigma_list,0.6) # step 1 切分数据集
    [sample_train,sample_test]=[lai_waiter_inc_sigma_list[:],lai_waiter_inc_sigma_list[:]]  # step 1 切分数据集
    
    logger.info("训练模型")
    logger.info("1.模型初始值")
    params_X0=[0.080,-6.985,1,1,0.771,-0.028] 
    logger.info("2. 训练模型")
    params_arr=train_WMCmodel(sample_train,params_X0,train_err_image_path)   
    logging.info("模型初值:\t{}".format(str(params_X0)))  
    logging.info("训练得到的模型系数:\t{}".format(str(params_arr)))
    logger.info("3. lai初始值     ")                                   
    lai_x0=np.mean(np.array([item[4] for item in sample_train])) 
    logger.info("4. 测试模型训练结果")
    err=test_WMCModel(sample_test,params_arr,lai_x0,test_err_image_path)    
    err_code=[i[3] for i in err]
    logger.info("4.1 剔除 误差大的点")
    err_max=10
    logger.info("阈值：{}".format(err_max))
    out_lines=np.where(np.abs(np.array(err)[:,1])>err_max)
    logger.info('误差点序：\t{}'.format(str(out_lines[0].tolist())))
    logger.info("4.2 保留误差点 较小值")
    in_lines=np.where(np.abs(np.array(err)[:,1])<=err_max)
    logger.info('保留点序：\t{}'.format(str(in_lines[0].tolist())))  
    logger.info("5. 模型评价")      
    err_arr=np.array(err)[in_lines,:][0,:,:]
    err_arr=np.array(err)                           
    model_evaluate(err_arr[:,0],err_arr[:,2],err_code,evalution_image_path)
    
    logger.info("生产影像") 
    result=process_tiff(sigma_tiff['data'].astype(np.float64),
                        incident_angle['data'].astype(np.float64),
                        soil_water_tiff['data'].astype(np.float64),
                        soil_simga_dbl['data'].astype(np.float64),
                        params_arr,lai_x0)
    np.save(os.path.join(result_dir_path,'lai_result.npy'),result)
    result.tofile("result.bin")
    logger.info("保存计算结果数据:\t{}".format(os.path.join(result_dir_path,'lai_result.npy')))
    imageHandler.write_img(os.path.join(result_dir_path,'lai_wmc.tiff'), sigma_tiff['proj'], sigma_tiff['geotrans'], result, no_data='null')
    logger.info("结果影像地址：\t{}".format(os.path.join(result_dir_path,'lai_wmc.tiff')))
    logger.info("程序运行结束")
